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Abstract 

Statistical data is often analyzed as a contingency table, sometimes 
with empty cells called zeros. Such sparse tables can be due to scarse 
observations classified in numerous categories, as for example in ge- 
netic association studies. Thus, classical independence tests involving 
Pearson's chi-square statistic Q or KuUback's minimum discrimination 
information statistic G cannot be applied because some of the expected 
frequencies are too small. More generally, we consider goodness of fit 
tests with composite hypotheses for sparse multinomial vectors and 
suggest simple corrections for Q and G that improve and generalize 
known procedures such as Ku's. We show that the corrected statistics 
share the same asymptotic distribution as the initial statistics. We 
produce Monte Carlo estimations for the type I and type II errors on a 
toy example. Finally, we apply the corrected statistics to independence 
tests on epidemiological and ecological data. 

1 Introduction and notations 

Physical, sociological or biological surveys often lead to data presented as 
contingency tables. These surveys aim at studying relationships such as 
total, mutual, partial or conditional independence between several characters 
in a population. Table notations are very useful to formulate the test and 
make the corresponding hypotheses explicit, but can be cumbersome when 
the number of characters exceeds three. For theoretical results, we therefore 
use vector notations instead, and we reformulate the independence test as a 
multinomial goodness of fit test. 



1.1 Goodness of fit tests 

Let p = {pi, . . . ,Pr) be a probability vector of dimension R where R ^ 2 
is the total number of cross-classifying categories. Let n be the sample 
size and x = (ni, . . . ,nfi) the vector of observed frequencies, realization of 
X = {Ni, . . . ,N{{) with distribution M{n;p), multinomial distribution with 
parameters n and p. We denote by p^ the probability vector under the null 
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hypothesis and consider the following test 



Ho- p = against Hi : p i= p^ ■ (1-1) 

Popular goodness of fit statistics include Pearson's chi-square statistic 
Q and KuUback's minimum discrimination information statistic G, defined 
in O [To]. For two probability distributions p and p' , these are written : 

Q,(p')=nX:^^^^, (1.2) 

and 

R I 

Gp{p') = 2nY,p'r^^-- (1-3) 

They belong to the power divergence statistics family {i?C^, A G M} defined 
by Read and Cressie in [5], respectively for A = 1 and A — )■ 0, where : 



9r7 ^ 
r=l 



/ \ A 

Pr 



For a review on goodness of fit testing methods and statistics, see [3]. 

The vector p^ is not always completely specified. We assume that p^ is 
a known function of an unknown parameter ^ of C R'' with s < R — 1, 
and denote p^ = p^{6) with 6 = {6i, ... ,6s)- We suppose that the functions 
9 I— >■ p^{6) we consider here are bijective and estimate p^{6) by p*^ = p^{9*) 
where 9* is the maximum likelihood estimator of 9. As the probability p 
is generally unknown, we also estimate the pr for r in {1, . . . , R} by their 
maximum likelihood estimators p* = rir/n. Throughout this paper, under- 
lying indexes n are omitted for simplicity of notation, and we consider the 
statistics Qp*a{p*) and Gp*o{p*). 

Read and Cressie show that under Birch's regularity conditions, see [2], 
the statistics RC^ are asymptotically equivalent and share a common chi- 
square limit distribution: 

Theorem 1. 

VAGM, hin Ch,{RC^o{p*)) = Xr-s-1- (1-5) 

This implies the following result : 
Corollary 1. 

lim Cho{Qp*o{p*)) = lim ^^Ho{Gp*o{p*)) = x%,^s^v (1-6) 

n— >-+oo n-^+oo 
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This asymptotic result is a consequence of the Central Limit Theorem. 
Classical empirical limitations include that the sample size n must be over 30 
and that all expected frequencies must be over 5. Usually Gpto{p*) is pre- 
ferred to Qp*o{p*) because it is less sensitive to small cell frequencies. Read 

and Cressie recommend the use of RC^io (p*) instead for n ^ 10 and a min- 
imum expected frequency over 1. Sometimes, even the lower bound 0.5 is 
accepted for the expected frequencies. A review on these conditions can be 
found in |4]. 

1.2 Sparse tables 

When there are too few subjects in the study or when the classifying cate- 
gories are too numerous, the table comprises one or several empty cells called 
random zeros. The table is then called sparse and it is likely that at least one 
cell has an expected frequency below 0.5. Random zeros would not appear 
if the sample was of sufficient size. Structural zeros, corresponding to cells 
with an expected probability of zero, are not considered here and should be 
suppressed. We therefore assume that the following condition is satisfied : 

p:V0, re {!,..., i?}. (1.7) 

Sparse tables can nonetheless be tested for independence, for example by 
regrouping cells so that the condition on the expected frequencies is satis- 
fied. However this procedure is not always data relevant. Fisher's exact test 
given in [6] applies without restrictions, except that it becomes numerically 
unmanageable when the table dimension grows. This leads to the use of 
Monte Carlo simulation methods as explained in |1]. 

The approach we propose here consists in correcting the historical statis- 
tics Qp*o{p*) and Gp*o{p*) according to the number of zero cells, by gener- 
alizing and improving a method designed by Ku in [9]. 

2 Corrections for Pearson's and Kullback's statis- 
tics 

Let C be the random variable giving the number of zeros in the vector 
or the contingency table, and c a realization of C with R — c ^ 1. For 
simplicity, we assume that ni = ?i2 = •••= nc = and that nj ^ 
1 for all j in {c -|- 1,...,R}. The maximum likelihood estimator p* = 
(0, . . . , 0, Uc+i/n, . . . , nji/n) underestimates the pi for i in {1, . . . , c} and 
overestimates the pj for j in {c + 1, . . . , R}. Its use when c 7^ thus has 
consequences on the statistics <5p*o(p*) and Gp*o{p*). 
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2.1 Ku's correction for one zero 

Ku argues in [9] that Gp.o(p*) tends to inflate with respect to Qp*o{p*) when 
C grows. He then proposes to subtract 1 from Gp*o{p*) for each zero, that 
is c in total. He proves the asymptotic equivalence between Gp*o{p*) and 
its corrected version only for c = 1. I will explain why his reasoning is 
inconsistent. First, he considers a new statistic which is not a member of 
the power divergence family. Moreover, he uses the straightforward Lemma 
[T]to deduce the approximation: 

(2-1) 

for a = nj./n, b = p*^ and = 0, despite the fact that a = and 1/a is not 
bounded. 

Lemma 1. For each a,b > such that a < 2b, b < 2a and the quantities 
1/a and 1/6 are bounded, we have : 

The expression on the left in (|2.ip is null whereas the one on the right 
is negative. The sum over r of the left-hand side gives Gp*o{p*) and we 
recognize Qp*o{p*) in the sum of the right-hand side. However, unlike Ku 
seems to think, zero and non-zero cells tend to compensate for each other 
and the approximation (|2.ip can not be extended to the corresponding sum. 
There is therefore no behavior of G and Q that we can deduce from this. 
Finally, he illustrates this asymptotic result on a small sample of size n = 10. 

We propose new corrections for both statistics, based on Ku's correction 
and a likelihood inequality that we present in the next subsection. 



2.2 Likelihood inequality 

Let us consider the following inequality coming from a likelihood reasoning : 
the sample vector we observe can be thought of more likely to happen than 
any other possible vector, since it is the one we actually observed. With 
exactly c zeros and R — c non-zero cells observed, so for all m ^ c and n'j 

such that n'j ^ nj, \/je{c+l,...,R} and n = ^'j + J2f=c+i ""i ■ 

P(iVi =0,...,N^ = 0,N,+i = n^+i, ...,NR = nR)^ 
F{Ni =n[,...,Nm = n'^, iV^+i = ... = N^ = {), N^+i = ...,Nr = n'j^). 

(2.2) 

We give in Proposition [1] a sufficient condition on the pr for ()2.2p to be 
satisfied, and prove this statement in Appendix |Aj 



4 



Proposition 1. The inequality (|2.2p is satisfied under the following assump- 
tion : 

Pi^^, Vie {l,...,c}, \/j e{c+l,...,R}. (2.3) 

n 

2.3 Corrections 

We propose an estimator p for p, different from the maximum likeliliood 
estimator, with the following form : 

(pi =a, Vi €{!,... ,c}, 

\p, =^-d, yj e{c+l,...,R}, ^ ■ ' 

where a = a^., b = bn and d = dn are random variables depending on n 
designed to compensate for the under- and overestimations due to p* . We 
thus take them positive with < 6 < 1, such that b inflates the modified 
maximum likelihood estimator nj/n^ and such that d controls the related 
rise. 

This new probability vector allows us to define corrected statistics Qp*o{p) 
and Gp*o{p), provided that the parameters a, b and d satisfy several con- 
ditions. Forcing the summation of the pr to 1 implies that (i? — c)d = 
ac + 7i^~^ — 1, and thus allows us to define Qp*o{p°'^) and Gp,o{p°'^) with j5"* 
such that : 



Definition 1. 



ac + n 



1-6 



Vi G {l,...,c}, 

1 (2.5) 
-, yj e{c + l,...,R}. 



R-c 

Note that for c = 0, we fix a = and 6=1. 

Before considering the other conditions on p""^, let us first give some 
notations. Let n be mmj^^^^i ji^{nj} and n be Ina,Xj^^c_^_l Jly{nj}. Let 
n stand for n — n{R — c) and n for n{R — c) — n. We dismiss the uniformly 
distributed case where : 

fi 

n = n = nj = — , e {c + 1, . . . , R}, (2.6) 

It — c 

thus guaranteeing that n and n are positive. Let 

/ Inj^/jR-l)) In(^) In(n-n) ^ 

6min = max 0, — . , . — , (2.7) 

\ [n[n) m(n) m(n) / 

a^in 6 = max 0, ^ . ' , 2.8 

V cn° I 



5 



and 



n'^ — n nf' — n 



a^Ub) = min 1, — ^, ^.(^(^ _ -) ^ ^) , (2-9) 



Proposition [T] is applied to p"^^. Together with inequahties < pf' < 1, for 
all r in {1, . . . , R} it is then equivalent to these conditions on a and b: 

bmin <b< 6max = 1, (2.10) 

and 

(b) <a< amax(fe). (2.11) 

We want to make a practical choice among possible values of p""^ ensuring 
us that it is as far from p* as possible. We therefore fix b in (|2.10p quite far 
from 1 as a convex combination of 6min and ^max with an empirical parameter 
h equal to 0.1. For this value of b we choose a near the upper limit of the 
interval in (|2.1ip . that is : 

b = ft.6max + (1 - h)bmm and a = amax(^) - e, (2.12) 

where e is a small constant designed to eliminate boundary effects. 

The final expressions we get for the corrected statistics Q"^^ = Qp*o{p°'^) 
and G"^ = Gp*o(p"^) of Q = Qp*o{p*) and G = Gp*o{p*) are: 

Q"'' = n2(i-^)Q-/(a,6), (2.13) 

with 



-a 



R — c ^ np*^ 

j=c+l J 

2 1 / ac + n — 1^ 



and 

where 



G^^ = n^^^G - g{a,b), (2.15) 



, ,ac + n^-^-l A, f ni(R-c)~n^(ac + n^-''-l)' 



(2.16) 
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2.4 Convergence 

In this paragraph, we show the convergence in distribution of Q"^ and G"''' 
to a chi-square distribution. Let us first study the parameter b. 

Property 1. Bound bmin is strictly less than 1. Moreover, if n = o[n) and 
n ^ n when n tends to +oo, then bmin o,nd b tend to 1. 

Proof. Inequahties n{R — c) < nR, n — n{R — c) < n and n — n<n show that 
all three components of the maximum defining 6min are strictly less than 1. 
Their order 1 developments as n tends to +cxd give their convergence to 1, 
with the quantities R — 1, R — c and R — c — 1 bounded. □ 

We also study the asymptotic behavior of C, for once denoted C„, and 
state the following lemma which proof appears in Appendix lAl 

Lemma 2. The number of zeros Cn converges almost surely to as n tends 
to +00. 

We recalled in section [1] the convergence of Pearson's and Kullback's 
statistics to a chi-square distribution. In Theorem [21 we state a similar 
property for the corrected statistics and G""^. 

Theorem 2. Under Birch's regularity conditions in J^, the estimator p""^ 
defined by (f23|) , (f2A0]) and (|2TT]) is such that : 

lim Cno{Q''')= lim Cn,{G'^') = xl-s-i- (2.17) 

Proof. We deduce from Lemma [2] the existence of a set 0,' of probability 1 
on which we can find a rank hq such that C„ = for n ^ hq. The variable a 
is then set equal to and the variable b equal to 1. Hence, the estimates p* 
and p""^ match, so that Q""^ = Q and G""^ = G. Theorem [1] then completes 
the proof. □ 

The final two sections are dedicated to proving the relevancy of our cor- 
rections through simulations and real data analyses. 



3 Simulations 

In this section, simulations confirm the necessity to correct not only G but 
also Q. We compute the statistics Q, G, RC^/^ = RG'^io{p*), Q"^ and G"'' 
on 1 000 vectors of length R = 100 of total frequency n = 400 for each of the 
four multinomial distributions defined in Table [T]by /i to /4. Let Er denote 
the expected frequencies for r in {1, ... , 100}. 

For each set of vectors sharing the same c, we give the quantile of order 
1 — a = 95% for the five statistics considered. Results are displayed in 
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Table 1: 


Multinomial probabilities /i to /4. 




|{r; < 0.5}| 


Probability 


fl 


20 


(0.0002, . . . , 0.0002, 0.01245, . . . , 0.01245) 

^ V ' ^ V ' 


/2 


50 


20 times 80 times 

(0.0002, . . . , 0.0002, 0.0198, . . . , 0.0198) 


/3 


70 


50 times 50 times 

(0.0002, . . . , 0.0002, 0.03286667, . . . , 0.03286667) 

^ V ' V ' 


h 


90 


70 times 30 times 

(0.0002, . . . , 0.0002, 0.0982, . . . , 0.0982) 

90 times 10 times 



figure [H as well as a line indicating the chi-square quantile of order 1 — a = 
95% with R—1 = 99 degrees of freedom Xo 95 99 = 123.22. Only the center of 
each graph should be considered because the quantiles for extreme numbers 
of zeros are computed on too few observations, sometimes only on one or 
two among the 1 000 simulations in total. 

Quantile values of Q tend to explode as c grows whereas G stays quite 
stable around Xo95 99- This behavior is the opposite of the one predicted 
by Ku. For /i, statistics Q"^ and G"'' lead to the rejection of the null 
hypothesis. For /2 to however, their quantiles lie below the critical line 
and Hq is accepted. We thus have compensated for the rise of Q, and both 
our corrected statistics are stable. 

This analysis is confirmed by the computation of empirical risks of type I 
for a = 0.01, 0.05 and 0.1 as showed in Table E] and by the power study 
below. Probabilities /i to are perturbed into /( to such that for all j 
in {1,2,3,4}: 

Vie{l,...,10}, /j(f) = /,(f) + 1/300, (3.1) 
Vie {11,..., 90}, f'j{i) = fj{i), (3.2) 
ViG{91,...,100}, /j(i) = /,(i) - 1/300. (3.3) 

Vectors are simulated with probabilities fi to f^, and goodness of fit for 
/( to is tested. The Tables |2] and |3] show that the empirical type I risks are 
lower for our corrections compared to the classical statistics when c is quite 
important, whereas the empirical power is much higher for our corrections 
when c is small. 

4 Applications 

We apply the total independence test using the corrected statistics Q°'^ and 
G"'^ to two datasets involving two-dimensional tables. For such tables, the 
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14 16 18 20 22 

c 








Figure 1: Quantiles of order 0.95 for Q, G, G"^ and RC^I^ as functions 
of c, under null multinomial probabilities /i to /4, for 1 000 samples of size 
n = 400 and R = 100 categories. The line represents the threshold 95 99- 
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Table 2: Empirical type I risks for Q, Q'^^ G, G"'' and RC'^/^, for 1 000 
samples of size n = 400 and R = 100 categories, at levels a = 0.01, 0.05, 
0.1, for vector probabilities /i to f^. 



a 


■Ho 


mode{c) 


Q 




G 


Qab 




0.01 


fi 


19 


0.031 


0.233 


0.003 


0.401 


0.003 




f2 


47 


0.030 


0.005 













h 


65 


0.027 
















h 


84 


0.031 














0.05 


h 


19 


0.044 


0.352 


0.010 


0.573 


0.017 






47 


0.024 


0.005 













h 


65 


0.069 











0.006 




h 


84 


0.052 














0.1 


h 


19 


0.130 


0.479 


0.033 


0.674 


0.039 






47 


0.072 


0.005 





0.005 







h 


65 


0.137 
















h 


84 


0.098 











0.006 



Table 3: Empirical powers for Q, G, G"^ and i^C^/^, for 1 000 samples 
of size n = 400 and R = 100 categories, at levels a = 0.05 for simulated 
vector probabilities /i to /4, and null vector probabilities /( to f'^. 



■Ho 


■Hi 


mode(c) 


Q 




G 


Qab 


RG^^ 


fi 


fi 


19 


0.233 


0.853 


0.322 


0.983 


0.157 


/2 


f'2 


47 


0.087 


0.026 


0.009 


0.061 


0.009 


/3 


f^ 


64 


0.229 


0.005 








0.027 


h 


fi 


84 


0.089 
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Table 4: Diplotype table for the association between TNFAIP3 and 

Systemic Sclerosis.. 



Status Hl/Hl H1/H2 H1/H3 H1/H4 H1/H5 H1/H6 

Sound 98 7 116 2 71 3 

Affected 91 9 104 3 70 12 

H2/H3 H2/H5 H2/H6 H3/H3 H3/H4 H3/H5 

Sound 4 2 34 1 42 

Affected 5 4 1 30 2 40 

H3/H6 H4/H5 H5/H5 H5/H6 

Sound 2 1 13 1 

Affected 7 1 13 5 



hypotheses are written : 

T-Lq: pij=pi+p+j against Hi: 3 {ii, ji) e I x J, pi^j^ pi^^p+j^, 

(4.1) 

where pi^ and p^j are the marginal distributions for the two characters 
featured in the table. To ensure the condition (|1.7p we remove the empty 
lines i of {1, ...,/} such that nj+ = 0, and the empty columns j of {1, . . . , J} 
such that n+j = 0. 

4.1 Multi-marker approach for Systemic Sclerosis 

Table |4]is the diplotype table obtained from an association study in Humans 
looking for an association between three genetic markers on the gene TN- 
FAIP3 and Systemic Sclerosis presented in |8]. Empty columns have been 
removed. A haplotype is the allelic distribution of markers on a chromosome, 
and a diplotype is the combination of both parental haplotypes. Diplotype 
tables, though more interesting than haplotypic tables because they take 
into account more information, are usually trickier to handle because they 
are sparse. Our corrected statistics can therefore be helpful in such situa- 
tions. 

Though nine haplotypes theoretically exist, only eight are observed, de- 
noted HI to H8, leading to 8^ = 64 diplotypes Hi/Hj. Two samples are 
compared, affected versus sound subjects, on which we test the indepen- 
dence between the diplotype configuration and the health status of n = 794 
individuals. 

The table is of dimension 2 x 16, that is i? = 32 categories with c = 1 zero 
and s = 16 parameters. There are exactly 16 expected frequencies below 5. 
The chi-square quantile Xq 95 15 = 24.99 is compared to the statistics : 

Q = 14.62, Q'^'' = 20.76, G = 15.82, G"'' = 28.43, RG"^/'-^ = 14.85. 
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Table 5: Contingency table for the joint study of trophic level and 
vegetable composition in rivers. 



Trophic level 


(0,0,0) 


(1,0,0) 


{r,P; 
(0,1,0) 


e) 

(0,0,1) 


(1,1,0) 


(0,1,1) 


Oligotrophic 








3 





3 


2 


Mesotrophic 


2 


1 





2 


1 





Eutrophic 


2 





3 


1 


1 






Only G leads to reject the null hypothesis of independence. This seems to 
be the right decision since it is confirmed by the single-markers approaches 
and haplotype tests in |5], all showing a significative association between the 
markers and the disease. 

4.2 Trophic level and vegetables in the rivers of the Petite 
Camargue Alsacienne 

The search for a link between the trophic level and the vegetable composition 
of some rivers of the Petite Camargue Alsacienne in North-East France leads 
to Table \5\ A river can be either oligotrophic, mesotrophic or eutrophic, if 
its nutritive content is respectively poor, intermediate or high. Uncommon 
vegetables are considered rare, exotic or polluo-tolerant. To each river, a 
triplet of binary characteristics (r, p, e) is assigned, indicating the presence 
(1) or the absence (0) of rare (r), exotic (e) and polluo-tolerant {p) species. 
The original ecological study can be found in |11) . We consider n = 21 
different rivers. Two empty columns were removed from the original table, 
leading to Table [5] of dimension 3x6, that is = 18 categories and c = 7 
zeros, with s = 17 parameters. Here are 3 expected frequencies below 0.5. 
Test statistics are compared to the chi-square quantile Xq 95 10 — 18.31: 

Q = 14.38, Q"'' = 20.68, G = 18.67, G"'' = 26.05, i^C^/^ = 14.84. 

Both corrected statistics Q""^ and G""^ as well as G lead to reject the null 
hypothesis, indicating an association between trophic level and vegetable 
composition. A thorough study of the table shows that rare species tend to 
settle preferentially in oligotrophic rivers. They are indeed better adapted 
to this kind of environment which tends to disappear from the rivers in 
the study. Moreover, polluo-tolerant species constitute the majority of the 
vegetables in eutrophic rivers. A eutrophic environment is competitive and 
these resistant species tend to get the best of it. 

4.3 Discussion 

We suggest to compute both Q"''^ and G"'^, and to reject the null hypothesis 
if at least one of them is larger than the chi-square threshold. 
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Our results tend to prove that this approach is relevant and our correc- 
tions efficient in sparse tables. They are all the more interesting for the fact 
that sparse tables are usually left aside because the hypotheses needed to 
apply classical chi-square tests are not satisfied. 



A Appendix section 

Proof of PropositionUl Assume that the ffi'st m of the c frequencies rij, 1 ^ 
i ^ c are modified. Let n'pj G {c + 1, ... , R} compensate for these changes. 
The likelihood inequality (|2.2p is then equivalent to : 

(nc+i-<+i) {nR-n' ) n[ n'^ 

Pc+1 Pr y,^^ ...P^. (A.l) 

(nc+i - <+i + 1)! (tir - ra^ + 1)! ^ n{\ n'^\ 

Let us show that (|2.3p implies (lA.ip . Applying (|2.3p to each element of the 

R m 

left hand side of (|A.ip with multiplicities such that : {uj —n'j) = n^, 

j=c+l i=l 

we get : 

x...x4^-^^p^...p< (A.2) 



As {rij - n'j + 1)! ^ n^''^""-J^ for all j in {c + 1, . . . , i?} and n'-l ^ 1 for all i 
in {!,..., m}, we deduce (lA.ip and equivalently (|2.2p . □ 

Proof of Lemma\^ Let us ffist show that: 

Ve > 0, lim P(C„ > e) = 0. (A.3) 

n— >+oo 

For each n we compute P(Cn = c) for c in {0, . . . , i? — 1}. Let be the 
probability under the null hypothesis = {p\, . . . ,p^^, p^_^_^, . . . ,p^). A 
subject belongs to one of the first c cells with probability q'^ = P? + • • • + Pc 
and to one of the R — c last cells with probability 1 — = p^_^i + ' ' ' + P/ji 
with g° in ]0, 1[. 

Let us consider now the binomial distribution B{n; q^) and write the 
probability P(C„, = c) of obtaining a table containing exactly c zeros placed 
anywhere : 

nCn = c) = (^^Y{Ni = ...=N, = 0, Nc+i^0,...,Nr^{(j^,4) 
= (A.5) 

Then : 

R-l 

¥{Cn = 0) = l-Y,nCn = c) = l-F{Cn>e), VeG]0,l[. (A.6) 



c=l 
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As ('l^) is bounded and (1 — g")" tends to 0, the probability P(C„, = c) also 
converges to for all 1 ^ c ^ -R — 1 when n tends to +oo, and so does the 
corresponding sum over c. 

We now use Borel-Cantelli's Lemma to conclude that C„ converges to 
almost surely. Indeed, for e > 0: 

R— 1 

5^ P(C„, >e)^Yl ^ 1) = E ( e ) i < (^-^^ 

n^l n^l c=l ^ ^ ^ 

□ 
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